function J = Jaco(x, y, z, a, b2, c, delta)
    % 差分法计算3D-SIMM系统雅可比矩阵
    x0 = [x; y; z];
    f0 = SIMM(x0(1), x0(2), x0(3), a, b2, c);
    J = zeros(3,3);
    for j = 1:3
        x_perturb = x0;
        x_perturb(j) = x_perturb(j) + delta;
        f1 = SIMM(x_perturb(1), x_perturb(2), x_perturb(3), a, b2, c);
        J(:,j) = (f1 - f0) / delta;
    end
end